Larval thermal performance curves
AS Huffmyer
2021
This script analyzes and plots data for Symbiontic Integration 2021 respirometry data for thermal performance curves.
Setup
Set up workspace, set options, and load required packages.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)## install packages if you dont already have them in your library
if ("tidyverse" %in% rownames(installed.packages()) == 'FALSE') install.packages('tidyverse')
if ("car" %in% rownames(installed.packages()) == 'FALSE') install.packages('car')
if ("scales" %in% rownames(installed.packages()) == 'FALSE') install.packages('scales')
if ("ggplot2" %in% rownames(installed.packages()) == 'FALSE') install.packages('ggplot2')
if ("nls.multstart" %in% rownames(installed.packages()) == 'FALSE') install.packages('nls.multstart')
if ("broom" %in% rownames(installed.packages()) == 'FALSE') install.packages('broom')
if ("rTPC" %in% rownames(installed.packages()) == 'FALSE') remotes::install_github("padpadpadpad/rTPC")
#load packages
library("ggplot2")
library("tidyverse")
library('car')
library('scales')
library('nls.multstart')
library('broom')
library('rTPC')
library('ggstatsplot')Data visualization
Load data from LoLinR.
PRdata<-read.csv("Pacu2021/Output/Respiration/TPC/oxygen_Resp_rates.csv") #load dataFormat data columns.
#remove all rows of wells that did not have samples or blanks
PRdata<-PRdata[!is.na(PRdata$Type),]
#format columns
PRdata$Temperature<-as.factor(PRdata$Temperature)
PRdata$Temperature<-as.numeric(as.character(PRdata$Temperature))
PRdata$Treatment<-as.factor(PRdata$Treatment)
PRdata$Plate<-as.factor(PRdata$Plate)
PRdata$Run<-as.factor(PRdata$Run)
PRdata <- PRdata %>%
group_by(Temperature,Treatment)%>%
mutate(abs=abs(R.nmol.org.min))%>%
summarise(rate=abs) %>%
mutate(temp=Temperature)
PRdata$group <- paste0(PRdata$Temperature,"_", PRdata$Treatment)View and remove outliers.
#PRdata <- PRdata %>% filter(R.nmol.org.min > -0.15) %>% filter(R.nmol.org.min < 0)
#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).outlier.plot
ggsave("Pacu2021/Figures/Respiration/TPC/outliersbygroup.pdf", outlier.plot, dpi=300, w=16, h=8, units="in")
q <- c(0.25, 0.75)
Quants <- PRdata %>%
group_by(Treatment, Temperature) %>%
summarize(quant25 = quantile(rate, probs = q[1]),
quant75 = quantile(rate, probs = q[2]),
IQRbyGroup=IQR(rate))
Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)
Quants$upper <- Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range
#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)
#remove outliers
PRdata <- PRdata %>%
filter(rate < upper) %>%
filter(rate > lower)
outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).outlier.removed.plot
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")Plotting
Plot data as a dot plot
r_plot<-PRdata %>%
ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) +
xlab("Temperature") +
scale_fill_manual(name="Treatment", values=c("blue","red"))+
scale_color_manual(name="Treatment", values=c("blue","red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter.pdf", r_plot, dpi=300, w=8, h=8, units="in")Filtering
PRdata <- PRdata %>%
filter(!Temperature==8) %>%
filter(!Temperature==18)%>%
filter(!Temperature==34)%>%
filter(!Temperature==36)%>%
filter(!Temperature==38)%>%
filter(!Temperature==40)
r_plot<-PRdata %>%
ggplot(., aes(x = Temperature, y = rate, colour=Treatment)) +
geom_point(aes(fill=Treatment, group=interaction(Temperature, Treatment)), pch = 21, size=2, alpha=0.5, position = position_jitterdodge(0.3)) +
#geom_smooth(method = "lm", se=FALSE, formula=y ~ poly(x, 3, raw=TRUE), size=3)+
xlab("Temperature") +
scale_fill_manual(name="Treatment", values=c("blue","red"))+
scale_color_manual(name="Treatment", values=c("blue","red"))+
ylab(expression(bold(paste("R (nmol ", O[2], " larva"^-1, "min"^-1, ")")))) +
theme_classic() +
theme(
legend.position="none",
axis.title=element_text(face="bold", size=16),
axis.text=element_text(size=12, color="black"),
legend.title=element_text(face="bold", size=14),
legend.text=element_text(size=12)
); r_plot
ggsave("Pacu2021/Figures/Respiration/TPC/Respiration_Scatter_Filtered.pdf", r_plot, dpi=300, w=8, h=8, units="in")View and remove outliers.
#identify outliers by Tremperature and Treatment groups
outlier.plot <- ggbetweenstats(PRdata,group, rate, outlier.tagging = TRUE)## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).outlier.plot
Quants <- PRdata %>%
group_by(Treatment, Temperature) %>%
summarize(quant25 = quantile(rate, probs = q[1]),
quant75 = quantile(rate, probs = q[2]),
IQRbyGroup=IQR(rate))
Quants$group <-paste0(Quants$Temperature,"_", Quants$Treatment)
Quants$upper <- Quants$quant75+1.5*Quants$IQRbyGroup # Upper Range
Quants$lower <- Quants$quant25-1.5*Quants$IQRbyGroup # Lower Range
#join outlier cutoffs with data
PRdata <- left_join(PRdata, Quants)
#remove outliers
PRdata <- PRdata %>%
filter(rate < upper) %>%
filter(rate > lower)
<<<<<<< HEAD
PRdata <- PRdata %>% filter(rate > 0.0015)
=======
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
PRdata <- PRdata %>% filter(rate < 0.14)
PRdata <- PRdata[-c(96,106),]
outlier.removed.plot <- ggbetweenstats(PRdata, group, rate, outlier.tagging = TRUE)## Warning: Number of labels is greater than default palette color count.
## Try using another color `palette` (and/or `package`).outlier.removed.plot
ggsave("Pacu2021/Figures/Respiration/TPC/Removed_outliersbygroup2.pdf", outlier.removed.plot, dpi=300, w=16, h=8, units="in")Analysis
TPC ftting
#Padifeld et al rTPC and nls.multstart: A new pipeline to fit thermal performance curves in r ### https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13585
# choose model
get_model_names()## [1] "beta_2012" "boatman_2017" "briere2_1999"
## [4] "delong_2017" "flinn_1991" "gaussian_1987"
## [7] "hinshelwood_1947" "joehnk_2008" "johnsonlewin_1946"
## [10] "kamykowski_1985" "lactin2_1995" "lrf_1991"
## [13] "modifiedgaussian_2006" "oneill_1972" "pawar_2018"
## [16] "quadratic_2008" "ratkowsky_1983" "rezende_2019"
## [19] "sharpeschoolfull_1981" "sharpeschoolhigh_1981" "sharpeschoollow_1981"
## [22] "spain_1982" "thomas_2012" "thomas_2017"
## [25] "weibull_1995"#sharpeschoolhigh_1981
# get start vals
start_vals <- get_start_vals(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
# get limits
low_lims <- get_lower_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
upper_lims <- get_upper_lims(PRdata$temp,PRdata$rate, model_name = 'sharpeschoolhigh_1981')
start_vals## r_tref e eh th
## 0.0456096 0.7027915 2.1208907 30.9411765## r_tref e eh th
## 0.04434088 0.76808284 2.42820170 30.94117647low_lims## r_tref e eh th
## 0 0 0 1upper_lims## r_tref e eh th
## 0.1232836 10.0000000 20.0000000 32.0000000# AMBIENT CURVE FIT
d.amb <- PRdata %>%
filter(Treatment=="Ambient")
# # get start vals
# amb.start_vals <- get_start_vals(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
#
# # get limits
# amb.low_lims <- get_lower_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.low_lims
# amb.upper_lims <- get_upper_lims(d.amb$temp,d.amb$rate, model_name = 'sharpeschoolhigh_1981')
# amb.start_vals
# amb.low_lims
# amb.upper_lims
amb.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.amb,
iter = 500,
start_lower = start_vals - 1,
start_upper = start_vals + 1,
lower = low_lims,
upper = upper_lims,
supp_errors = 'Y')
amb.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
<<<<<<< HEAD
## 0.06369 0.77099 8.81372 31.66205
## residual sum-of-squares: 0.03744
##
## Number of iterations to convergence: 25
=======
## 0.06287 0.81631 10.05527 31.49888
## residual sum-of-squares: 0.03564
##
## Number of iterations to convergence: 21
>>>>>>> 2935501d0ff9c1aaecd004f8a449a38a9358f695
## Achieved convergence tolerance: 1.49e-08# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.fit, newdata = amb_new_data)
amb.TCP.res <- calc_params(amb.fit) %>%
mutate_all(round, 2) # round
amb.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.55 0.81 34.96 0.7 2.31 2.5 5.41 34.15
## breadth skewness
## 1 4.99 -1.61## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.07 29.58 2.34 34.35 0.77 2.9 2.71 4.77 32.01
## breadth skewness
## 1 4.58 -2.14BOOTSTRAPPING AMBIENT CURVE FIT
Resampling the original data with replacement
# refit model using nlsLM
amb.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.amb,
start = coef(amb.fit),
lower = low_lims,
upper = upper_lims,
weights = rep(1, times = nrow(d.amb)))
# bootstrap using case resampling
amb.boot1 <- Boot(amb.fit_nlsLM, method = 'case')
# look at the data
head(amb.boot1$t)## r_tref e eh th
<<<<<<< HEAD
## [1,] 0.05622679 1.0115372 3.274674 32.00000
## [2,] 0.06621550 0.6028136 12.494659 31.84514
## [3,] 0.11625218 1.5155506 4.932409 29.19690
## [4,] 0.12328360 1.6623770 3.525770 28.08559
## [5,] 0.12328360 1.5686641 2.828922 27.31338
## [6,] 0.10076807 1.2078562 3.945774 29.13973# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot 
# HIGH CURVE FIT
d.high <- PRdata %>%
filter(Treatment=="High")
# choose model
mod <- 'sharpschoolhigh_1981'
# # get start vals
# high.start_vals <- get_start_vals(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
#
# # get limits
# high.low_lims <- get_lower_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.low_lims
# high.upper_lims <- get_upper_lims(d.high$temp,d.high$rate, model_name = 'sharpeschoolhigh_1981')
# high.start_vals
# high.low_lims
# high.upper_lims
high.fit <- nls_multstart(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.high,
iter = 500,
start_lower = start_vals - 1,
start_upper = start_vals + 1,
lower = low_lims,
upper = upper_lims,
supp_errors = 'Y')
high.fit## Nonlinear regression model
## model: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 28)
## data: data
## r_tref e eh th
## 0.1233 1.4038 2.4282 26.4288
## residual sum-of-squares: 0.03235
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.fit, newdata = high_new_data)
high.TCP.res <- calc_params(high.fit) %>%
mutate_all(round, 2) # round for easy viewing
high.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.44 7.45 57.05 0.53 0.47 2.01 29.62 49.6
## breadth skewness
## 1 8.87 0.06BOOTSTRAPPING HIGH CURVE FIT
Resampling the original data with replacement
# refit model using nlsLM
high.fit_nlsLM <- minpack.lm::nlsLM(rate~sharpeschoolhigh_1981(temp = temp, r_tref,e,eh,th, tref = 28),
data = d.high,
start = coef(high.fit),
lower = low_lims,
upper = upper_lims,
weights = rep(1, times = nrow(d.high)))
# bootstrap using case resampling
<<<<<<< HEAD
high.boot1 <- Boot(high.fit_nlsLM, method = 'case')
# look at the data
head(high.boot1$t)## r_tref e eh th
## [1,] 0.1232836 1.0875773 1.3645045 24.82696
## [2,] 0.1232836 1.0663405 1.7193422 25.18904
## [3,] 0.1232836 1.5667784 3.4832088 27.14641
## [4,] 0.1232836 1.2804056 2.9616468 26.35147
## [5,] 0.1232836 1.5242106 3.1892471 27.27148
## [6,] 0.1232836 0.0798746 0.1515831 1.00000##
## Number of bootstraps was 997 out of 999 attempted# look at the data
head(high.boot1$t)## r_tref e eh th
## [1,] 0.09466776 1.1017541 3.245382 28.74252
## [2,] 0.12328360 1.5033158 2.816571 26.55534
## [3,] 0.12328360 1.3592890 2.435695 26.91997
## [4,] 0.12328360 1.6925379 2.945170 26.25989
## [5,] 0.07216997 0.8133763 4.153968 29.78028
## [6,] 0.04920720 0.3859032 4.080503 32.00000# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = sharpeschoolhigh_1981(temp, r_tref, e, eh, th, tref = 28))
# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
high.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), high.preds, col = 'red') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot
# plot data and model fit
TPC.plot <- ggplot() +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Metabolic rate')
TPC.plot 
ggsave("Pacu2021/Figures/Respiration/TPC/TPC.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_SharpSchool.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")#confidence intervals of TPC parameters
broom::tidy(amb.fit_nlsLM)## # A tibble: 4 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
<<<<<<< HEAD
## 1 r_tref 0.0637 0.0110 5.78 5.07e- 7
## 2 e 0.771 0.412 1.87 6.76e- 2
## 3 eh 8.81 7.46 1.18 2.43e- 1
## 4 th 31.7 0.850 37.2 1.34e-37broom::tidy(high.fit_nlsLM)## # A tibble: 4 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 r_tref 0.123 0.497 0.248 0.805
## 2 e 1.40 3.76 0.373 0.710
## 3 eh 2.43 1.14 2.12 0.0384
## 4 th 26.4 20.7 1.28 0.208#AMBIENT
amb.extra_params <- calc_params(amb.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
amb.ci_extra_params <- Boot(amb.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')##
<<<<<<< HEAD
## Number of bootstraps was 83 out of 200 attemptedamb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"
#HIGH
high.extra_params <- calc_params(high.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
high.ci_extra_params <- Boot(high.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')##
<<<<<<< HEAD
## Number of bootstraps was 139 out of 200 attemptedhigh.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"
#Join Ambient and High estimates and CIs
All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>%
mutate_if(is.numeric, round, 2)
#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("blue","red"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
facet_wrap(~param, scales = 'free_y') +
scale_x_discrete('')
estimate.plots
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_SharpSchool.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")#Gaussian
#AMBIENT
amb.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.amb,
iter = c(4,4,4),
start_lower = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') - 1,
start_upper = get_start_vals(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987') + 1,
lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
supp_errors = 'Y',
convergence_count = FALSE)
amb.gaus.fit ## Nonlinear regression model
## model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
## data: data
## rmax topt a
## 0.06347 27.83799 4.38800
## residual sum-of-squares: 0.03694
##
## Number of iterations to convergence: 11
## Achieved convergence tolerance: 1.49e-08# predict new data
amb_new_data <- data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), 0.5))
amb.preds <- augment(amb.gaus.fit, newdata = amb_new_data)
amb.TCP.res <- calc_params(amb.gaus.fit) %>%
mutate_all(round, 2) # round for easy viewing
amb.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84 17.1 38.58 1.37 0.73 NA 10.74 21.48
## breadth skewness
## 1 5.87 0.65amb.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.amb,
start = coef(amb.gaus.fit),
lower = get_lower_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.amb$temp, d.amb$rate, model_name = 'gaussian_1987'),
weights = rep(1, times = nrow(d.amb)))
# bootstrap using case resampling
amb.boot1 <- Boot(amb.gaus.fit_nlsLM, method = 'case')
# look at the data
head(amb.boot1$t)## rmax topt a
## [1,] 0.06197087 27.13551 4.854462
## [2,] 0.05900555 27.60181 4.673603
## [3,] 0.05973728 28.41771 4.843638
## [4,] 0.06388289 27.78184 4.564628
## [5,] 0.05356135 27.06235 5.075532
## [6,] 0.07453052 28.19593 3.481052# create predictions of each bootstrapped model
amb.boot1_preds <- amb.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.amb$temp), max(d.amb$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))
# calculate bootstrapped confidence intervals
amb.boot1_conf_preds <- group_by(amb.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
amb.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
amb.CI.plot 
# HIGH
high.gaus.fit <- nls_multstart(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.high,
iter = c(4,4,4),
start_lower = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') - 1,
start_upper = get_start_vals(d.high$temp, d.high$rate, model_name = 'gaussian_1987') + 1,
lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
supp_errors = 'Y',
convergence_count = FALSE)
high.gaus.fit ## Nonlinear regression model
## model: rate ~ gaussian_1987(temp = temp, rmax, topt, a)
## data: data
## rmax topt a
## 0.04746 27.07125 5.77867
## residual sum-of-squares: 0.0274
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 1.49e-08# predict new data
high_new_data <- data.frame(temp = seq(min(d.high$temp), max(d.high$temp), 0.5))
high.preds <- augment(high.gaus.fit, newdata = high_new_data)
high.TCP.res <- calc_params(high.gaus.fit) %>%
mutate_all(round, 2) # round for easy viewing
high.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01 14.15 28.29
## breadth skewness
## 1 7.72 -0.26high.gaus.fit_nlsLM <- minpack.lm::nlsLM(rate~gaussian_1987(temp = temp, rmax, topt, a),
data = d.high,
start = coef(high.gaus.fit),
lower = get_lower_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
upper = get_upper_lims(d.high$temp, d.high$rate, model_name = 'gaussian_1987'),
weights = rep(1, times = nrow(d.high)))
# bootstrap using case resampling
high.boot1 <- Boot(high.gaus.fit_nlsLM, method = 'case')
# look at the data
head(high.boot1$t)## rmax topt a
## [1,] 0.05927776 26.86467 4.041476
## [2,] 0.05146475 26.79550 4.483799
## [3,] 0.04765019 30.66779 8.565635
## [4,] 0.04808687 28.83889 8.571844
## [5,] 0.04426714 25.74916 8.505839
## [6,] 0.04689024 28.19428 6.521274# create predictions of each bootstrapped model
high.boot1_preds <- high.boot1$t %>%
as.data.frame() %>%
drop_na() %>%
mutate(iter = 1:n()) %>%
group_by_all() %>%
do(data.frame(temp = seq(min(d.high$temp), max(d.high$temp), length.out = 100))) %>%
ungroup() %>%
mutate(pred = gaussian_1987(temp = temp, rmax, topt, a))
# calculate bootstrapped confidence intervals
high.boot1_conf_preds <- group_by(high.boot1_preds, temp) %>%
summarise(conf_lower = quantile(pred, 0.025),
conf_upper = quantile(pred, 0.975)) %>%
ungroup()
# plot bootstrapped CIs
high.CI.plot <- ggplot() +
geom_line(aes(temp, .fitted), high.preds, col = 'red') +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Respiration Rate (nmol O2/larva/min')
high.CI.plot 
# plot data and model fit
TPC.plot <- ggplot() +
geom_point(aes(temp, rate), d.amb, size = 2, alpha = 0.5,col = 'blue') +
geom_point(aes(temp, rate), d.high, size = 2, alpha = 0.5,col = 'red') +
geom_line(aes(temp, .fitted), amb.preds, col = 'blue', size=2) +
geom_line(aes(temp, .fitted), high.preds, col = 'red', size=2) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), high.boot1_conf_preds, fill = 'red', alpha = 0.3) +
geom_ribbon(aes(temp, ymin = conf_lower, ymax = conf_upper), amb.boot1_conf_preds, fill = 'blue', alpha = 0.3) +
theme_bw(base_size = 12) +
labs(x = 'Temperature (ºC)',
y = 'Metabolic rate')
TPC.plot 
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_Gaussian.pdf", TPC.plot, dpi=300, w=8, h=8, units="in")
amb.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.06 27.84 17.1 38.58 1.37 0.73 NA 10.74 21.48
## breadth skewness
## 1 5.87 0.65high.TCP.res ## rmax topt ctmin ctmax e eh q10 thermal_safety_margin thermal_tolerance
## 1 0.05 27.07 12.93 41.22 0.53 0.79 2.01 14.15 28.29
## breadth skewness
## 1 7.72 -0.26#confidence intervals of Gaussian TPC parameters
broom::tidy(amb.gaus.fit_nlsLM)## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rmax 0.0635 0.00604 10.5 2.38e-14
## 2 topt 27.8 0.570 48.9 1.67e-44
## 3 a 4.39 0.778 5.64 7.48e- 7broom::tidy(high.gaus.fit_nlsLM)## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 rmax 0.0475 0.00492 9.65 2.91e-13
## 2 topt 27.1 0.847 32.0 2.56e-36
## 3 a 5.78 1.60 3.61 6.71e- 4#AMBIENT
amb.extra_params <- calc_params(amb.gaus.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
amb.ci_extra_params <- Boot(amb.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(amb.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')##
## Number of bootstraps was 109 out of 200 attemptedamb.ci_extra_params <- left_join(amb.ci_extra_params, amb.extra_params)
amb.ci_extra_params$Treatment <- "Ambient"
# Joining, by = "param"
#HIGH
high.extra_params <- calc_params(high.gaus.fit_nlsLM) %>%
pivot_longer(everything(), names_to = 'param', values_to = 'estimate')
high.ci_extra_params <- Boot(high.gaus.fit_nlsLM, f = function(x){unlist(calc_params(x))}, labels = names(calc_params(high.gaus.fit_nlsLM)), R = 200, method = 'case') %>%
confint(., method = 'bca') %>%
as.data.frame() %>%
rename(conf_lower = 1, conf_upper = 2) %>%
rownames_to_column(., var = 'param') %>%
mutate(method = 'case bootstrap')##
## Number of bootstraps was 159 out of 200 attemptedhigh.ci_extra_params <- left_join(high.ci_extra_params, high.extra_params)
high.ci_extra_params$Treatment <- "High"
#Join Ambient and High estimates and CIs
All_params <- rbind(amb.ci_extra_params, high.ci_extra_params)
All_params <- All_params %>%
mutate_if(is.numeric, round, 2)
#Plot
estimate.plots <- ggplot(All_params, aes(Treatment, estimate, color=Treatment)) +
geom_point(size = 2) +
scale_color_manual(name="Treatment", values=c("blue","red"))+
geom_linerange(aes(ymin = conf_lower, ymax = conf_upper)) +
theme_bw() +
facet_wrap(~param, scales = 'free_y') +
scale_x_discrete('')
estimate.plots
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")
ggsave("Pacu2021/Figures/Respiration/TPC/TPC_estimates_Gaussian.pdf", estimate.plots, dpi=300, w=8, h=8, units="in")